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I. INTRODUCTION 



One of the long standing paradoxical problems in theoretical physics lies in understanding 
how a macroscopic system reaches equilibrium departing from the reversible microscopical 
laws of nature. The simplest systems for which the origin of irreversibility can be studied on 
a microscopic basis are the linear ones 1T|,0]. In this work we investigate this behavior for a 
Brownian particle in a quantum-mechanical heat bath composed of a finite number of small 
oscillators. The key point of this analysis is the consideration of the linear coupling between 
the Brownian particle and the bath, which allows us to reduce the Hamiltonian to a set of 
uncoupled oscillators. This model is recurrently studied in the literature from many different 
approaches, such as the Langevin equation, the master equation, and the exact solution for 
the evolution operator, using different techniques. The importance of this model lies in 
its broad applicability in many fields of physics: condensed matter, statistical mechanics, 
quantum optics, quantum electrodynamics, quantum measurement theory, scattering and 
decay theory, etc. In the majority of previous works statistical fluctuations, dissipation, 
and equilibrium tendency in linear quantum-mechanical systems are shown to result from a 
projection of the total quantum system onto a restricted subspace. A macroscopic equation 
is obtained corresponding to a reduced description of the system. The restriction to the 
Markovian approximation and weak-coupling limit are usually undertaken as well as the 
analysis is carried out in the limit of a dense bath. Only few works are devoted to study 
the behavior of the system for a finite (discrete spectrum) bath P,^. In this case it is 
not possible to prove convergence to an equilibrium state in the limit t ^ oo because of 
the existence of Poincare recurrences, which become extremely infrequent for large systems. 
This is the reason why previous works prefer to eliminate them by passing to the limit of 
an infinite heat bath. However, quantitative estimates of the Poincare recurrence times for 
finite systems can be made in a way compatible with a dissipative behavior. In this work 
we reinforce this fact. We study the time evolution of a finite system departing from the 
exact solution of the eigenvalue problem without appealing to approximations, assuming a 
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factorizable initial condition where the bath has reached a unique thermal equilibrium state 
(passive reservoir). No use is made of coarse graining, finite memory assumptions, randomly 
varying Hamiltonians or nonlinear modifications to the Schrodinger evolution. We perform 
our calculations for an arbitrary spectral density and temperature. 

The work is organized as follows. In Sec. II we introduce the model and its exact solution. 
We also show a criterion to disregard a broad class of coupling functions. The exact solution 
is used in Sec. Ill in order to study the time evolution of the relevant variables of our 
problem, i.e. the mean occupation number (or energy) and mean position of the subsystem, 
the behavior of bath variables, etc. Sec. IV is devoted to derive an exact generalized form 
of the quantum Langevin equation, with time-dependent coefficients. Numerical results and 
their analysis are presented in Sec. V. We show that the mean position operator performs 
damping oscillations correlated with the mean energy of the subsystem, which decays in 
time towards a state of equilibrium with the bath. After reaching this state revivals occurs 
periodically in the subsystem. The bath remains almost unaltered in thermal equilibrium as 
a consequence of its passivity and robustness. In Sec. VI we take the limit of a dense bath and 
obtain, by means of an analytic continuation method, a complex frequency Zq — D,+5D,+ir/2 
into which the unperturbed real frequency ft of the oscillator is shifted by the heat bath. 
For a long period of time the exponential decay law dominates the evolution. In this period 
the standard form of the Langevin equation is derived, where the mean displacement of 
the Brownian particle undergoes a slowly damped harmonic oscillation corresponding to the 
complex frequency Zq. Deviations from the exponential decay law are also discussed. 

II. BROWNIAN MOTION: THE MODEL AND ITS EXACT SOLUTION 

Let us consider a harmonic oscillator interacting with a bath modeled by a set of harmonic 
oscillators. The Hamiltonian of the system is 




(1) 
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where Hj represents the interaction. Capital and lower-case letters stand for subsystem 
and bath variables respectively. In our case Hj only involves a linear coupling between the 
Brownian particle and the bath, i.e. 

where all c„ are real and small coupling constants. We define, as usual, creation and anni- 
hilation operators {% = 1) 



(3) 



2 ' "Y 2m„un^"'' 

that satisfy the canonical commutation relations 



h W 



(4) 



the other commutators vanish. In terms of these operators the Hamiltonian reads 



H = n [B^B + \ + Y.^n [hihr. + i,)+T.9n [Bhl + B^) , (5) 



n=l ^ ^' n=\ 

where gn = Cn/\/ Mflmn^Jn- The linear interaction allows us to find normal modes of H 
in an exact way. This kind of coupling is known in the literature as the rotating wave 
approximation |^ (in general only the coupling between coordinates is taken into account). 

In this model the interaction term preserves the total number of quanta. In fact, defining 
the number of quanta operators 

Nn. = B^B, 



the total number of quanta given by 



Nn = bibn, 
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N 
n=l 

is a constant of motion, due to dNx/dt = i [Hj, N^] = 0. Then we can resolve the Hamil- 
tonian into sectors of definite number of quanta. For A^^ = 1 (one-particle sector) and 
calling 

\n) =B^\o) = I0...0), 

(6) 

W =&J.|o) = |o)®|o..._L...o), 

we obtain 

N N 

Hi = n \n) {n\ + J2iUn \iOn) + E ^« (1^) Kl + 1^") (^1) + (7) 



n=l n=l 



TV 



where C = ^ + J2 This is the discrete version of the Friedrichs model [Q. 

n=l 

Let us find now the normal modes of H, i.e. the new set of uncoupled harmonic oscillators 
with normal frequencies a,y (Greek subscripts run from to A^, while Arabic ones run from 
1 to A^). Hence we write H as 



N 



H = Y, a,4c, + C, (8) 

i/=0 

where the new creation operators, Cy, are related to the old ones by means of a unitary 
(canonical) transformation 

N 

= ^^B + C^^nbn, (9) 
n=l 

which preserves the canonical commutation relations 

V- (10) 



Coefficients 
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are the matrix elements of the unitary change of basis, from {!«;/)} to , \ujn)}, where 

= 4 |o) 

is an eigenvector of Hi 

N 

Hi = ^a,y \a^) {a„\ + C. 
The canonical commutators (p!0|) impose the following condition for and (^^n 

N 

+ (P^nKn = ^M., (11) 
n=l 

which in the one-particle sector is a consequence of the orthogonality among eigenvectors of 
Hi, 

By taking into account transformation (P) in the Heisenberg equation of motion for Cy, 

= [cy,H] = a^Cy, (12) 



dt 

and using 



N 

i^ = [B,H] = nB+ Egnbn, 

n=l 



= [bn,H] =U}nbn + gnB, 



we obtain the following system of linear equations 



N 
n=l 



From the second of these equations we can obtain 0y„ as 

(Pun ) 
Otu — LOn 
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(13) 



(14) 



which is vahd only if ^ uJn-, Vz/, n. Replacing it into the first equation of (|T3D we have 

= 0. (15) 



Since we are looking for non-trivial solutions the expression between brackets must be iden- 
tically zero. Then we have an equation for the normal frequencies of the new set of harmonic 
oscillators 



N 



a,. 



9n 



n=l - 

The procedure developed above is equivalent to solve the eigenvalue problem E.\ |q;^) 
ay \ay) for the matrix that represents the Friedrichs Hamiltonian in the basis , Icj^)} 



(16) 



/ 



^ g\ ■ ■ ■ qn 
■ 







9n 



(17) 



In fact (plTf ) is the complex conjugated matrix of the matrix form of Eqs. ([T3|) because gn is 
real. From Eq. ([TTl) for fi = u [normalization of the eigenvectors of (0)] we obtain 

1 



TV 



(18) 



1+ y ( 

n=l 

which can be completely determined if we know the set of eigenvalues a^. The normal 
frequencies can be obtained numerically or by analytic perturbative methods (in some 
special cases can even be obtained exactly), so we assume that they are well known. In Fig. 
1 we show where these values are located. 
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FIG. 1. Location of the normal frequencies. 



The normal frequencies correspond to the intersection of the straight line a,^ — Q and 
the summation of the hyperboles g'^/ {a^ — Un) ■ From this picture we see that the normal 
frequencies always lie between consecutive frequencies of the unperturbed Hamiltonian, ex- 
cept for the two extremum values which lie outside the interval delimited by cui and cun- 
{q:u}u=o,...,n never coincide with {wn}n=i,...,N, and are very close to each cUn near the extrema 
and move away from uun in the centrum. In order to guarantee the positivity of the Hamil- 

N 2 

tonian as the lowest frequency approaches zero and the convergence of the series E — , 

n=l "'"-'^ri 

it is required that 

N 2 N 2 

V ^ r <n-u;i + d and V % < con + S - Q, (19) 

where S is an infinitesimal parameter (e.g. the distance between contiguous unperturbed 
frequencies), and also that g{u!) behaves smoothly around cui and cun, having a small value 
at these points. Then the coupling privileges the interaction with the subsystem oscillator 
of frequency fl. These conditions express the fact that the interaction is small. However, for 



small N, the interaction must not decrease very fast from the centrum since if such were the 
case the subsystem oscillator would be coupled to few bath oscillators and then the bath 
would not be effective. These problems disappear approaching to the continuum (Sec. VI). 



III. TIME EVOLUTION 



We are looking for the way in which the relevant variables of the system evolve in time. 
That is to know the time evolution of B and 6„ and, from them, all other related dynamical 
variables. Coming back to Eq. (|12D we can integrate it and obtain the temporal evolution 

of C,y 



Using the closure relation 



c,{t) = e-*"''*c,(0). (20) 



N 



u=0 

we can obtain the following identities: 



N 

E . = 1, 

N 

E $:0.n = 0, (21) 

N 
u=0 



Eqs. ( pTD allow us to perform the inverse of the transformation 



N 



B = Y: Kcu, (22) 



N 

bn = T.€nCu. (23) 



u=0 

We are interested in knowing the explicit form of B{t) and &„(t). Then, from Eqs. ( pOD 
(0), (15), and d^) we have 
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N N 

B{t) = E ^^e-^^-'^c^O) = E $*e-*"'^* 

u=0 u=0 



N N 

bn{t) = E 0:„e— ^*c,(0) = E 0:„e— ^* 



N 



$,5(0) + E 0.n&„(O) 



n=l 



TV 



(24) 



$.5(0) + E 0.™&m(o) 

m=l 



We have expressed the time evolution of the unperturbed annihilation operators only in 
terms of their initial values. By considering Eq. (|14|) we write Eq. (p^) as 



AT 



B{t) = E I^^.Te"^""* 

u=0 



N 



(25) 



N 



hJt) = E \^uV^^e 



1^=0 



TV 



^(0)+ E ;r^&r.(0) 



in=l 



and a similar set of equations for the Hermitian conjugate operators. Eqs. ( p5|) are deter- 
mined uniquely by knowing the eigenvalues a^. 

These are the exact solutions of our problem. From them we can obtain all the relevant 
information under consideration. For example, one of the relevant variables is the position 



of the subsystem oscillator, X{t) = 

N 



B^{t) + B{t)] , given by 



X{t) = Y,\'^uf{ [cos (a,t) X(0) + sin (a^t) P(0) 



u=0 



(26) 



N 

E 



where P = Similarly, from P = iJ ^ 



, , , sin iayt) 

y/rrinUJri COS [a^t) Xn{0) H Pn{0) 



B\t)-B{t) 



we have 



_ AT 

P(t) = ^ I [_ sin (a,t) X(0) + cos (a^t) P(0) 



(27) 



+ 



N 

E 



COS (ttjyt) 

Vm„t^„ sm [ayt) g„(0) H , . Pn(0) 



Another interesting magnitude is the occupation number of the oscillator representing 
the Brownian particle (subsystem dynamics). We have 
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N 



N 

B'B) (0) + E [Bh\) (0) 



AT 



+ E 



9n 



n=l - 



AT 

(5t6„)(o)+ 5: 



n=l - 



(28) 



(0) 



Let us first consider the case in which the set of harmonic oscillators modeling the bath is 
in thermal equilibrium inside a big reservoir and the Brownian oscillator is isolated from the 
rest. At t = we extract the bath from the reservoir at temperature T, and put it in contact 
with the subsystem oscillator, in such a way that the bath becomes a thermal reservoir for 
the Brownian particle. In this situation the initial state of the total system is represented 
by a time-independent density matrix which is a direct (tensorial) product of the matrices 
representing the isolated harmonic oscillator Pb(0) and the environment degrees of freedom 
Pfe(O) in thermal equilibrium at temperature T 



p(0) = pb(0) 



(29) 



where Hi, = \ + |) is the bath Hamiltonian and tr^ is the partial trace over the 

n=l ^ ' 

reservoir. If it is the case we have no correlations among the initial states of subsystem and 
reservoir. To obtain the time evolution of {(b^ B) {t)) = tr {p(0) [b^b) (t)}, we need to 
specify the initial values of Eq. (pSj) in the state p of Eq. ([29|) . They are given by 

[b^b) (0)] 



{{Bbi) (0)) = = ((i?t6„) (0) 



(30) 



where k is the initial number of quanta in the subsystem oscillator. Therefore the subsystem 
dynamics is given by 



N 

{B^B)it)) = 2|$j'|$,|'cos[(a^-«,)t] 



N 



N 



K 



+ E 



9n 



(31) 



TV 



In a similar way we obtain the mean value of the number of quanta operator for the 
n— oscillator of the bath, i.e. 



N 



K 



+ E 



9m 



N 



\ (a^ - u;m)(ai/ - (^m) e^'^'" - 1 

(32) 



u=0 



{ay - iOn) 



N 



m=l 



The expressions obtained can be formally rewritten as 

N 

{Na{t)) = Pnn{t) (iVn(O)) + ^ P^nit) (iVn(O)) , 

n=l 

(33) 

N 

(7V„(i)) = P„n(i) (A^f^(O)) + Yl Pnm{t) (iVm(O)) , 

m=l 

where (Nq) — (-B^-B^ and (A^„) = {b\firi^ . In a forthcoming paper we will show that this is a 
general result of this kind of models, which allows us to derive the Pauh master equation. It 
can be proved that P^n and Pa„ are respectively, the transition probability of the one-particle 
state |Q) remaining unchanged (survival probability) and the transition probability from the 
state \u}n) to the state They represent the probability that at time t the contribution 
to the oscillator occupation number comes from itself and from the bath, respectively. Pno. 
and Pnm are the probability that the n— bath occupation number has contribution from 
the oscillator and from the bath, respectively. These probabilities satisfy the normalization 
condition 

TV TV 

Pq,q. + X/ -^^^ ~ ^' -^"^ X/ -^""^ ~ ^' 

n=l m=l 

and are explicitly given by 
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n 



AHt 



N 



N 



Pnn(t) — Pnuit) 



-iHt 



fi.i'—O 



OJr. 



u=0 



N 



2^1$ 1^1$ 1^ 9lcos[{a^- a^)t] 



fijiy — Q 



N 



(34) 



-iHt 



LOr. 



N 



2 E I^mI I'^.I 



aWm cos [{a^ - au) t] 



l-L,l' — 



N 
i/=0 



9n9n 



1 2 



We can see from the second equation of (^) that, although there is no interaction term 
in the Hamiltonian among the bath oscillators themselves, the time evolution for a bath 
oscillator has contributions coming from the whole bath. This fact was noticed in Ref. 0] 
[cf. Eqs. (2.2e) and (2.2f)]. 

The decomposition made in Eq. ( P^ ) is useful for studying the different contributions to 
the time evolution of the mean number operators (Sec. V). 

In the limit of low temperatures T ^ 0, (A^n(O)) = (e^""" - -> 0, and then {Nn{t)) = 
Pnnit) {Nfi{0)) = K,Pnfi{t). In the case k = 1 it is the survival probability (the probability 
of no decay of the state \^)), 



{Nam 



T=0 



n 



iHt 



n 



for {Nnm = 1. 



(35) 



In Sec. VI for a dense bath we show that the asymptotic behavior of this probability obeys 
a power-law decay. This is a well known fact in decay theory of unstable quantum systems 
and was reported as an anomaly in statistical treatments of quantum open systems (see, 
e.g., Ref. [§). 

In the next section we derive an exact equation of motion of the mean value of the 
position operator X (a generalized form of the Langevin equation). 
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IV. LANGEVIN EQUATION 



Let us consider Eqs. (^) and (p7|) for the bath in thermal equihbrium at the initial 
time. In this case let {|A^n)} be a basis of eigenvectors of A^„. Taking mean values in the 
state (E^ we have 



{bnm = tr 



N 



PB(0)exp[-/3a;„(iV„ + l/2)] U exp[-f3uj,{N, + 1/2)] 

tr^e-^^"} 



There is a vanishing factor 

ENjNnlbniO) exp[-Pun{Nn + l/2)]\Nn) , siucc (iV„ | (0) |iV„) = 0. Then (6„(0)) = 0. 
Similarly (bi{0)) = 0. So we have (g„(0)) = = (p„(0)) . Thus 



N 



(X(t)) =E|$.r cosKt)(X(0))+sinKt)(P(0) 



N 



(36) 



P{i)) =El'^-r -sinKt)(X(0)) + cosKt)(P(0) 



v=0 



The initial and instantaneous variables are related by a generalized sum of rotations 



Pit) 



N 



u=0 



I 



V 



cos{aut) sin(ajyt) 



a transformation which can be summarized as 



m 



a{t) h{t) 
-bit) a{t) 



\ 



( 



V 



(^(0)) 
P(0) 



/ 
V 



(^(0)) 

m 



(37) 



(38) 



where a{t) = J2 \^u\ cos{ai^t) and b(t) = \^u\ sin (0;,^^) . 

We can invert the matrix of Eq. (^) to obtain (X(0)) and ^P(O)^ as functions of {X{t)) 
and (P{t)) : 



(39) 



^ (X(0)) 1 


1 


^a(t) 


-bit) ^ 


^ (xit)) ^ 


V (^(0)) ) 


~ A 


^b{t) 


a{t) ^ 
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where A(t) = a^(t) + By taking time derivatives in Eq. ( pSD and replacing the initial 
mean values by those of Eq. ( p^ ) we have 



V 



Pit) 



J 



1 

A 



a a+ bo b cl— a b 



d b— b a d a+ b b 



V 



(Xit)) 

Pit) 



(40) 



Similarly, from the second derivative of (|38D we have 

d a+ b b b a— b 
d b— b a d a+ b b 









A 

V 



V 



(xit)) 
Pit) 



(41) 



Finally eliminating (^Pit)j from Eq. ( ^OD we obtain a generalized form of the Langevin 
equation with time-dependent coefficients, 



X it)) + n^t) {Xit)) + Tit) (x it)) = 



(42) 



where 



n\t) 



d b - b a 



b d —a b 



(43) 



a b —b d ' ' a b —b d 

The standard stochastic force /stoch does not appear in Eq. (|4^) since it is included in the 
terms containing the operators g„(0) and Pn(0), which were eliminated by taking the mean 
values in a thermal equilibrium initial state of the bath. That is (/stoch (^)) = 0. Eq. ( ^21) 
contrasts with the equivalent, but non-local in time, standard integro-differential form of 
the equation of motion of (X) . Eq. (^2]) is actually a rather complicated expression since 
the coefficients are not easy of evaluating. In Sec. VI we estimate them in the continuous 
limit. 



V. EXAMPLES AND RESULTS 

Let us describe the model we used for obtaining the numerical results of this section. 
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A. Choice of parameters 



The model described in Sees. II and III consists of three main ingredients: the subsystem 
and the bath, the interaction, and the initial conditions. We have considered a subsystem 
represented by a harmonic oscillator with natural frequency Q and mass M (heavy Brownian 
particle), a bath of small oscillators with frequencies c<j„ varying in a range between u^in 
and Climax, a small linear coupling between system and bath, gn = Ac^, where A = 
and c'„ = c„ (m^u;^)"^''^ , and the whole composed system prepared in such a way that at 
t = 0~ there is no correlation between subsystem and bath. The bath is in equilibrium with 
an external heat source at temperature T = , where kB is the Boltzmann constant, 

and at t = 0"*" the bath is extracted from the thermal source, put in contact with the 
subsystem and the total system is left isolated. Each parameter mentioned above defines a 
typical time scale. These are: the scale associated with the natural frequency of the isolated 
subsystem, the scale defined by the lowest frequency of the bath, uj^l^, related to the 
reaction of the system when the interaction is switched on; the decay time F"^ [see Eq. (69) 
in Sec. VI] in which the subsystem dissipates its energy into the reservoir, and which is 
related to the squared of the perturbation parameter A; the memory time related with the 
highest frequency present in the bath, uj^l^; the time scale P associated with thermal effects 
(relative to quantum ones); the Poincare recurrence time given by the minimal difference 
between contiguous normal frequencies, specifically tp min(a^+i-a ) t^^^ ^'i- (^^) below]; 
and two time scales related with quantum deviations form the exponential decay law, a very 
short time tz (Zeno period, which is responsible for no decaying of the subsystem under a 
continuous succession of measurements and occurs because of the temporal derivative of the 
survival amplitude vanishes at t = 0) and a very long one tx (Khalfin period of power series 
tails, which is a consequence of the lower bound of the energy) (see Sec. VI). In order to 
have a manifestation of these time scales the parameters and variables must be chosen with 
certain criterion. An important condition we must take into account and which is frequently 
overlooked in the literature is condition (|19|). For example, in the case of a semi-infinite 
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frequency spectrum, oj G (0, oo), the often used ohmic spectral density, g^ioj) ~ does 
not satisfy condition (|19|), since it has a logarithmic type divergency. As we want to obtain 
a pictorial image of the temporal evolution of the main magnitudes of Sec. Ill, then we 
specify the parameters appearing in these magnitudes [e.g. Eqs. p^), (p^), (pT]), and (p2D] 
as follows: 

VL = l, 

(44) 

K=l. 

This choice of (3 fixes the thermal time scale to the same value of the VL~^ scale, and then 
purely quantum-mechanical and thermal effects are comparable, so we are far of the classical 
limit hVt <^ ksT. The choice of k facilitates the comprehension of the one-particle sector 
several times studied in decay theory. 

For the sake of simplicity we consider for the variables the case in which the bath fre- 
quencies are equidistant around the frequency Q, i.e. 

ujn = ^ + A(n - ^ , n = l,...,N, (45) 



where A is the spacing between contiguous frequencies of the bath, A = Un+i —ujn, being the 
band width cjjv — uji = A{N — 2), and the number of small oscillators is an odd integer. 
The coupling function is given by a Lorentzian-like function 

9n = -—- pT2- (46) 

This function is plotted in Fig. 2. D is related to the coupling strength A and is taken equal 
to A, for reasons which will become clear below. We fix a = in order to have half of 

the maximum value of gn at the extrema. Finally we take the band width equal to 0.018, 
for all A^, which is the value that allows us to compare our numerical results with those 
obtained by Gruver et al. [Q (who solved a set of coupled differential equations coming from 
a maximum entropy principle approach) in the case A^ -|- 1 = 32. In their work the choice 
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of the coupling function is different and also different the criterion to increase the number 
of small oscillators. While we maintain fixed the band width, they maintain fixed the value 
of the spacing A. As one of the purposes of this work is to study the way of reaching the 
continuous limit (see Sec. VI), we must have A —>■ 0, N —>■ oo, and AN =const. 
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FIG. 2. The coupling function. 



1.010 



Let us see the criterion used for fixing D. It is a consequence of conditions (|T^). Taking 
5 = A we have 



N 

E 



9l 



< ujn + A - il. 



^-^UJn + A - Un 

From Eq. (^) we have cu^ + A — Vl = ^^^^^^ + A — ujn = A{N + 1 — n), and for a large 
bath we can approximate 



16[{N/2f+{N/2-nf 



8^2 ^ 



=1 (iV + 1 - n) 



r. Then it must be satisfied 



< 1. 



1 2 



The summation can be bounded by times its maximum value, which is reached for n = N. 
Thus 
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D''N'', 1 



2A^ 



< 1, 



which imphes 



D < V2A. 



(47) 



So we find an upper boundary for D. In general we choose D = A. 

The + 1 normal frequencies a^, are obtained from a matrix-diagonalization routine. 

B. Numerical results 



In Figs. 3 to 7 we plot (Nq) vs. Qt. We show that (A^f^) decays in time to an asymptotic 
value, when the subsystem reaches equilibrium with the bath, given by (e — 1)"^ ^ 0.582 
[see Eq. (77) in Sec. VI, recall l3Q = 1]. After a long period a revival appears reaching 
again a similar value to the initial condition. Since for positive values of t the arguments 



of the cosines in Eq. (^) never are exactly in phase, the revival does not fully reconstruct 
the initial condition, so that the peak is smaller than the initial one and slightly broadens. 
Then, in the continuous limit, the time of revival goes to infinity and the peak gets out 
of sight among thermal fiuctuations. As shown in Fig. 9, this revival is periodic in time. 
It corresponds to the Poincare recurrence time and is given by the inverse of the smallest 
difference of normal frequencies, since Eq. ( ^1]) is a quasi-Fourier series because of the 
quasi-equidistance of the normal frequencies. Specifically 

2n 



(48) 



mm [^a^+i - a„) 

where obviously min (a^ — Oj,) = min (ttj^+i — a^) ■ Table I shows this time for different 
values of A^ + 1. 

Tab. I. Poincare recurrence time. 



A^ + 1 


10 


32 


100 


500 


tp 


3370 


11190 


37311 


177994 
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Fig. 3 shows that for a small number of bath oscillators, it does not result effective 
and loses the necessary robustness to break the natural oscillations of the subsystem. For 
+ 1 = 10 (Fig. 4) the energy lost of the subsystem oscillator is close to a dissipative 
behavior and the bath begins to be effective. It is surprising that for so few bath oscillators 
the subsystem already dissipates. With increasing A^, fluctuations get smaller and, since 
the spacing between frequencies decreases, tp grows. In Fig. 8 we draw {Nq) for different 
values oi N + 1 vs. a, re-scaled time with respect to tp {Q — 1). It shows the tendency to 
an exponential decay when the model approaches to the continuum. Nevertheless in the 
continuous limit the exponential decay law is not exact as we show in Sec. VI. 




0.5 I ' 1 ' 1 ' 1 ' 1 ' 1 ' 1 0.5 I ' 1 ' 1 ' 1 ' 1 ' 1 ' 1 

1000 2000 3000 4000 5000 6000 1000 2000 3000 4000 5000 6000 

fit fit 



FIG. 3. + 1 = 6. FIG. 4. + 1 = 10. 
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FIG. 5. N + 1 = 32. 



FIG. 6. iV + 1 = 100. 



50000 100000 150000 200000 250000 

fit 




FIG. 7. + 1 = 500. 



FIG. 8. Exponential approximation. 



Figs. 9 and 10 show the behavior of {Nq) for very long times in order to see the periodicity 
of this magnitude and the recurrence oi tp. In Fig. 9 we see that the height of the peaks 
monotonously decreases and afterwards it begins to oscillate. We do not have an explanation 



21 



of this fact. For even longer times we see that there does not exist a definite pattern repeating 
itself (Fig. 10, notice that in this picture only the envelopement of the peaks is plotted). 
Fig. 11 shows the form of a peak (A^ + 1 = 32) which is non-symmetrical. The growing side 
of the peak is steeper than the subsequent quasi-exponential decay. In Fig. 12 we choose 
D = 20A. In this case conditions (|TU|) are not satisfied and a non-dissipative behavior occurs 
with very quick oscillations. For Figs. 13 and 14 we have selected the value D = 2A for 
+ 1 = 32 and 100 respectively. This value satisfies conditions (0) [remember that the 
bound ( ^7|) is excessive]. However this worsening of D is reflected in the fact that {N^) 
presents more fluctuations. 




iJ A I 



2000000 4000000 6000000 8000000 10000000 



FIG. 9. Behvior of the peaks of {Nn). 



FIG. 10. Peaks for long times. 
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FIG. 11. Form of a peak. 



FIG. 12. N + 1 = 32, D = 2QA. 
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FIG. 13. + 1 = 32, D = 2A. 



FIG. 14. AT + 1 = 100, D = 2A. 



In Figs. 15 to 18 we plot {X) superposed to {Nq) for times around the second peak 
centered in tp. We see how the behavior of (X) and (Nq) are correlated and how after the 
revival the subsystem oscillator is damped [see Eq. (74) in Sec. VI]. In Fig. 16 the choice 
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of D — 2A shows again the growth of the fluctuations. 




Fig. 19 contains the behavior of Pqq and P^n for + 1 = 32. We see that the bath 
contribution remains almost constant except at times around tp in which the survival prob- 
ability has a maximum. In Fig. 20 we compare the behavior of (Nn) for a value near to the 
central frequency Q and for a value far of it. We see that while (A^2) smoothly fluctuates 
around a constant value (^e^^^ — 1^ , (Ni^) is sensible to what happens with {Nq) and then 
it is displaced with respect to (^e^^^^ — 1^ . This is an indication that the transference of 
energy from the subsystem to the bath is more effective for frequencies near to il. In Fig. 21 
we confirm the hypothesis that the central oscillators are those which receive the energy of 
the Brownian particle, since going to the continuous limit the distribution of P^n approaches 
to a delta function. The asymptotic value of (A^n) is given by 

N -, 

(iV^(oo)) = E^f^n(oo)-^^^, (49) 

n=l ^ ^ 

where Pq„(oo) = Yl, \\^v\ ^ ) = (^") • ■'■^ continuum we show, in Sec. VI, that 
Pq.u]{,oo) = 6{iJ — Q) up to the first order in a re-scaled parameter r = XH ^ (which is 
known as X^t approximation). Fig. 21 thus plots On vs. 




FIG. 19. Survival probability and bath contribution. FIG. 20. Bath population. 
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FIG. 21. Energy transfer to the central bath oscillator. 

We end this section by obtaining the results of Ref. for + 1 = 32. In Fig. 22 we 
show (Nq) vs. Qt (cf. Fig. 1 (b) of Ref. [Q). In Fig. 23 we see the damping oscillations, in 
Fig. 24 we compare the survival probability with the bath contribution, and in Fig. 25 we 
show the behavior of different values of (A^„) . 
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FIG. 22. (iVa) vs. nt. 



FIG. 23. {X) vs. nt. 
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FIG. 24. Survival probability and bath contribution. FIG. 25. (A^„) us. fit. 
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VI. ASYMPTOTIC LIMIT AND CONTINUOUS BATH 



In Sec. V we have studied a particular model for a finite but increasing number of 
bath oscillators with the aim of approaching to the continuous limit. In this section we 
analytically take this limit for frequencies spanning a segment into the positive real axis. 
At first glance the continuous limit of Eq. (|31|) can be made by changing summations for 
integrals. However it cannot be straightforwardly taken because of the appearance of a 
continuous set of divergencies along the integration domain. We must give a criterion to 
avoid these singularities, which is related to the choice of the boundary conditions. Following 
UUersma's pioneering work we propose a way, based on an analytic continuation method, 
to do that. 

We can define the function Rd{z) (the reduced resolvent operator in the energy repre- 
sentation, where d stands for discrete case) of the complex variable z departing from Eq. 



(|T^) as 



A'' 2 

R-,\z) = z-n-Y.^^. (50) 



where the normal frequencies a^, are given by the simple poles of Rd : R^^i^a^) = 0. Eq. (p!8[) 
can be rewritten in terms of Rd as 



= ^-^^ • (51) 



Eq. (|5lD allows us to write the first equation of Eq. ( p5|) as 



B{t) = E 



i/=0 



N 



B{0) + E ^^&n(0) 



„=1 - UJr, 



In order to perform the continuous limit let us consider the following identities: 



(52) 



y -T^, E KM = — f dz^^ E ^^&n(0), (54) 
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where C is a counterclockwise contour in the z-plane that encircles the + 1 singularities 
of Rd in the positive real axis (see Fig. 26). 



Im(z) 




^ — X — X — X — X- 




Re(z) 



FIG. 26. Counterclockwise contour C. 



We make use of the residues theorem: §q ^^^^ — ^tt^ Y.k=i ^^^[ 



2ni E 



Pjzk) 

fc=l Q'(^fc)' 



where are the simple zeros of Q{z). Eq. (^) is a direct consequence of the residues theorem 



and Eq. 



follows from the same theorem and from the fact that J2 



n n-l 



9nhn{^) = 0, 



since i?^^(co'„) diverges for all n. 

When the bath frequencies form a dense set, the normal frequencies are also dense. This 
limit of a continuous bath is valid for times that satisfy 

t <^ min (tti^+i — a^,)'^ . 



In this approximation R^^{z) goes to 

R-\z) = z-n- 



duj 



Z — ijJ 



(55) 



where cUmin and c^max are the extrema of the continuous set (lower and upper cutoff respec- 
tively) and avoid infrared and ultraviolet divergencies respectively. 
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g'^{uj) is defined by 



The function R~^{z) lias a cut along (cUmin, i^max), corresponding to tlie continuous spectrum 
of normal frequencies. In order to ensure that the equation R~^{z) = has no real roots it is 
necessary that Q G (cUmin, i^max), which together with conditions ([19|) adequately generalized 
to this case, provide a necessary and sufficient criterion of feasibility for dissipation in linear 



models. Generalized conditions ([T9|) are given by 



duj < ^7 — CUminj / duo < i^ra&yi — ^, (56) 

assuming that the integrals are well defined. The cut of the function R^^{z) can be reached 
from above and below the positive real axis, giving the limiting values ± ie) = a — 

^ ~ a-!li±ie ' where a G (tUmin, t^max)- Then, contraction of the contour C in Eq. (|53D 

yields 

^1 1 



1 ^ ~izt 1 

' -dz= / dae-'''' 



2TriJc R-^(z) 2-Ki 



R^^(a-ie) R~^(a + ie) 



(57) 



On the other hand, by taking into account the well known identity between distributions 

^ PV-^in6{x), (58) 



X ±ie X 

where PV stands for the principal value and 6 is the Dirac delta distribution, we have 

R~\a + ie) - R-\a - ie) = du ^ ^ . - du ^ ^ ' , = 2mg^{a). (59) 

Jajmin a — uj — le Jujn^in a — UJ + 16 

Then using ( p9D Eq. ( |57D is reduced to 

1 / e~*^* /■^m^"' fi'^(a) 



-dz= da ^e-"^\ (60) 

2mJcR~\z) \R-\a + it)\^ ^ ' 

where we consider that R^^[a — it) = R^^*{a + ie). Contracting the contour C in Eq. 
and taking the continuous limit we have 
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/ dz 



-izt 



2m Jc R-Hz) 



z — uj 2m 



dae 



-iat 



dujgioj) 



X 



1 



R~^[a — ie) {a — oj — ie) R~^{a + ie) {a — lu + ie) 
and after performing a straightforward calculation we obtain 



i> dz ^ 



-izt 



2m Jc R~Hz) 



du^b^O) 

Z ~ UJ 



da 



R~^{a-ie) 



e-*"*6„(0) 



+ / da duj 2 

-mm "''^min \R~'^{a + ie)\ [a — u + ie) 



-iat 



bUO), 



where we have used identity (0) for a; = a — a; to have the (5-function expressed as 

1 1 



{a — LU — ie) {a — LU + ie) 
Then B{t) can be written in the continuous limit as 



2TTi5{a — to). 



Bit) 



dae 



-iat 



i$,i'5(o) + <i>: 



dlU(l)a{^^)bu{0) 



where 



(61) 



(62) 



R-^(a + ie) 



6„(u;) = 5{a - Lu) + 



9{a)9i^) 



R ^(a + ie) (a — UJ -\- ie) 



(63) 



(64) 



Compare Eqs. ( |63D and (|6^ with those obtained in Ref. [Eqs. (4.8a) and (4.8b); in this 
work cjmin = and cUmax = oo], which correspond to the Lippmann-Schwinger coefficients of 
the eigenvectors of the continuous generalization of the one-particle Hamiltonian (|^) p|,p!0 



We have obtained the continuous generalization of B(t). By an straightforward but sim- 
ilar calculation we can obtain the continuous version of bn{t). 

From Eq. (|6^) we can give an approximate expression of the Langevin equation. By 
taking mean values in a thermal initial state for the bath we have 
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{B{t))= e-*"*(B(0)). (65) 

Performing an analytical continuation to the complex plane, one can extract the contribution 
of the poles of |$q,| . Then, let us analyze the analytic structure of as a complex 

function. The function R~^{z) has no zeros in the complex plane. It can be easily seen by 
considering z = a + ih, 



a + ib — UJ Ji^u.in (o- — + fe^ 

+ib{l+ duo- ^ ^ ' 



[a - a;)2 + 62 y 

where the imaginary part is equal to zero only if 6 = 0. Thus only real zeroes can exist, but 
in this case R~^{z) is only well defined by its limiting values R~^{a±ie). The discontinuity of 
R~^{z) in (ciJniin, cUmax) IS giveu by Eq. (|59D. _R~^(a±2e) has no zeroes because the imaginary 
part is not null, since g^i^a) 7^ for a G (cUmin, i^max) • Then R^^{z), which is analytic in 
the complex plane except for the cut discontinuity along (t^min, f^max) , has no zeroes in its 
definition range. In order to define an analytic function in all the complex plane we must 
analytically extend R~^{z) into the second Riemann sheet. To see this let us consider the 
function f{z) = {z — uj)^^ with its limiting values {a — 00 ± ie)~^. We can consider f{z) as a 
multivalued function or take the limiting values as defining two different functions (complex 
distributions). Then we define 



UJ 

1 



, for Im^; ^0, 



— =F 2712(5(2; — cu), for Im2;^0. 



for a E (t^mm,tUmax) , (^6) 



> 

f±{z) is analytic in all the complex plane. Let us check this fact for f+{z). Approaching the 
real axis from above and below we have 

1 1 

/+(« + ie) — /+(a — ie) = — h 2ni6{a — u), 

a — UJ + It a — UJ — It 

but that is exactly zero due to identity (^B|). Putting the functions f±{z) into R~^{z) we get 



two analytic continuations, Rj^i^z) (see, e.g., Ref. ||rT|). We can analytically continue the 
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function g'^{a) across (cUmin, ^^max) • We call G{z) this extension, i.e. there is an open region 
A of the complex plane containing (cUmin? "-^max) and a meromorphic function G : A C such 
that g'^{a) = G{a) for a G (ci^min, (^max) • For notational convenience we write g^{z) = G{z) 
also for non-real z. So, we are now ready to use the Cauchy theorem and change the contour 
of integration in Eq. by a contour E in the lower complex plane as shown in Fig. 27, 
which leaves the singularities of G{z) outside, namely 



k 

lm(z) 




FIG. 27. Contour E. 

The last part of definition ( |BBD can be considered as the values of f{z) into the second 
sheet. The relevant fact is that R'^^{z) has now a complex zero formally given by the zero 
of 

z-n- dcu^^-^ + 2inG{z) = 0. 

J'^aiin Z — UJ 

This zero can be estimated up to the second order as 
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Zq = fl 



duj ^ ^ ^ ' + 2mG{n) = n + 5n-i-. 

il — CO — ie 2 



That is, Q has two corrections, a real one, which provides the frequency shift, given by 



6n = pv 



duj- 



Q — u 

and a negative imaginary part, which represents the frequency width, given by 

r = 27^/(^]). 



(68) 



(69) 



For the model we have developed in Sec. V this second order correction should vanish due 
to the symmetrical distribution of both the bath frequencies and the interaction g{uj). 

Coming back to our original purpose let us evaluate the integral of Eq. (^7|) for the case 
in which only a simple zero zq is present. By taking the contour S lying below zo and using 
the residues theorem we have 



(Bit)) , 

R+) 1^0 

where the tilde over Rz^iz) stands for 



-izQt 



-iat 



R+\a)RZ\a) 



(70) 



being 6{a — Zq) the complex extension of the Dirac delta, defined by 

nax 

daf{a)S{a - zq) = f{zo), if zq e IntCg, 



which actually means dz' ^^^^ = f{zo), for Cq a contour encircling zq into the second 

sheet as depicted by Fig. 28. 
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Re(z) 



FIG. 28. Contour Cq; dotted line encircles zq into the second sheet. 



The second term inside the brackets of Eq. ([70|) is called the background. It is responsible 
for deviations from the exponential law and its contribution is relevant only for either very 
short or very long times. In the regime where the exponential decay dominates we can 
neglect the background. Then we obtain the approximate expression 



{B{t)) 



which satisfies the following differential equation 



{Bm 



(71) 



B m + z', {B{t)) = 0. 



(72) 



Eq. (|72|) has the form of a harmonic oscillator equation, however zq is now complex. Keeping 
in mind that zq = fl + 6Q — i^, and neglecting the terms in F^, it is derived 



B (t)) + + 6nf {B{t)) + t(b (t)) = 



(73) 



35 



which corresponds to the equation of a damped harmonic oscillator. Making the same for 
the creation operator and since X = l/\/2MQ [b + B"^^ we finally reach the standard form 
of the Langevin equation (in mean values): 



X (t)) + (fi + 5^Y {X{t)) + r (X it)) = 0. (74) 

The first two terms on the left hand side represents the Hamiltonian evolution of an oscil- 
lator with the renormalized frequency Vt + while the third one represents the 'friction' 
(dissipative) part of the damped linear oscillator with the damping factor F. This F is the 
one we could have obtained if we had performed the same analytical continuation made in 
this section and if we had neglected the background in Eq. (43). 



As we have pointed out in Eq. (|35|) the survival probability of the state \VL) is given 
by the square modulus of the integral of Eq. (|65D. Let us study its short and long time 
behavior. From general grounds, if we have a Hamiltonian given hy H = Hq + V, where 
Hq is an unperturbed Hamiltonian and V a small perturbation, and {iV^n)} is a set of 
eigenvectors of Hq, Hq \ipn) = En lipn) , the survival probability of the state {ipk) is given 
by Pk{t) = {ipkl e~*^Vfc) • At very short times the exponential can be approximated by 
its first two terms in the Taylor expansion, e"*''^* ^ 1 — iHt. Then the survival probability 
behaves as Pk{t) ~ 1 + O(t^), which does not correspond to an exponential decay behavior 
like 1 — Ft. Thus, at very short times, we have a non-exponential behavior, as it was shown 
is Sec. V, which is known as Zeno's period [P]. For very long times we also have a non- 
exponential contribution to the survival probability. As t goes to infinity the integral of Eq. 
(|65|) (survival amplitude) goes to zero as a consequence of the Riemann-Lebesgue theorem. 
Then the behavior of the survival probability depends on the small-frequency behavior of 
g'^ia). For small frequencies i?"^(u;min + «e) = t^min - ^ - J^"!"" (iu; ^^^J,. ~ uj^in - ^, 

mill '^min I ^ 

since the integral is bounded by condition (^6]). The behavior of g{a) is model dependent. 



We consider Ullersma's spectral strength of the kind g(a) = ^= — , where ci and C2 are 
constants, so the small-frequency behavior is given by g{a) ^ (ci/c2)a. Therefore we have, 
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for ^— < t < 



da|<l>„|'e-*°*^ - / da( ^-f: ~ / ^^"a e" 

n 



min 



Calling at = u we obtain for the survival amplitude 



•Tnax n ■ . o /'^maxt n • o / o • ■ • \ I^maxt 



Since the squared modulus behaves as t ^ [("^max ~ ^min) + 4] , the survival probability 
gives a power law decay 

Pnn{t) ~ AH-^ as t oo, (75) 



where A = y ci^max ~ '^min is a measure of an upper cutoff frequency. This deviation from 
the exponential decay law given by a power series tail is known from Khalfin's original work 

HI. 



Let us now study the asymptotic behavior of {NQ{t)) in the case of a dense bath, using 
the X^t approximation. Going back to Eq. (^) and its complex conjugate we write 

{Nn{t)) = da tiaV("-")* (iVf,(0)) 

•J UJinin '^'^min 

+ / da rf«V(°-° ^VaC' / c/a;0:H0,,H(iV^(O)), 



min " '-*-'min 



where we have considered an uncorrelated initial state between subsystem and bath, with 
the bath in thermal equilibrium, i.e. (A'^(O)) = (^e^^ — . From the Riemann-Lebesgue 
theorem all the oscillating terms vanish as t — > oo, so the only survival term is that with 
the product of deltas contained in 0's coefficients [Eq. (^4D]; namely 

/t'Jmax n'^(a) 1 

d^ 'l\^e^- 1 - ^^^^ 

We now perform the limit A — > 0. In such a limit R^^{a) ~ a — Q ± ie and taking into 
account Eq. (|59|) we obtain 
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(iVn(oo)) = / da— 



(77) 



_a — — ie a — ri + ie_ 

since the expression between square brackets is 2ni6{a — Q). We then see that the subsystem 
oscillator reaches thermal equilibrium with the bath. In Sec. V we have studied the behavior 
of the corresponding asymptotic form of (Nq) for a finite size bath [Eq. (^Of )]. In that case 
we have seen that, when approaching to the continuous limit, J2^=o {^ a^-i^ ) behaves as a 
delta function. It is important to remark that Eq. (|77|) shows that the transfer of energy 
mainly occurs between two oscillators in resonance. In Fig. 29 we have a diagrammatic 
representation of the time evolution of {NQ{t)) . First, there is a deviation from the expo- 
nential decay law due to the Zeno period, after that the exponential decay dominates for 
long time the evolution until the Khalfin power series tail, finishing in the asymptotic value 
of thermal equilibrium given by Eq. ( [77D . Note that the continuous limit carries Zeno time 
to zero and Khalfin time to infinity but nevertheless the exponential decay law is not valid 
at all. Numerical estimates will be given elsewhere. 
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VII. CONCLUDING REMARKS 



This work provides an exhaustive analysis of the most popular model of Brownian motion 
in a way which has not been deeply explored in the literature on the subject till now. 
The exact solution of the eigenvalue problem allows us to study the time behavior of the 
magnitudes of interest without resorting to approximations. No doubts about numerical 
errors can arise, since the diagonalization method used, which has a powerful speed of 
calculus, does not have recursive increasing deviations. Moreover the continuous limit is 
performed in an analytic manner, obtaining the standard results found in the literature. 
Figures of Sec. V clearly show all the properties expected for a dissipative system (a damped 
oscillator in this case), with estimates of the Poincare recurrence time, fluctuations and 
equilibrium. The Poincare period arises as an exact time of revival and not as an statistical 
property of the ensemble. The validity of the exponential decay law is also enlightened. In a 
forthcoming paper we will study other statistical properties of the model, such as correlation 
functions, and we will consider other spectral densities and different ways to distribute the 
unperturbed frequencies of the bath oscillators. 

On the other hand, the problem of irreversibility can be traced in the following way. Even 
for a finite system one can objectively 'see' an irreversible evolution, which only depends on 
the system and not on the ability of the observer. This irreversibility is not a consequence of a 
coarse-grained distribution or due to approximations. Nevertheless the time evolution is not 
strictly irreversible but it is practically irreversible for our scale of observation. Moreover, 
if we consider a real system with a large number of degrees of freedom (e.g. Avogadro's 
number) , it is easy to convince ourselves that we will see a time asymmetrical evolution for 
the Brownian particle, since the Poincare time becomes larger than the age of the Universe. 



39 



ACKNOWLEDGMENTS 



We are grateful to the organizers of the First International Colloquium on 'Actual Prob- 
lems in Quantum Mechanics, Cosmology, and the Primordial Universe' and the 'Foyer 
d'Humanisme' for their warm hospitality in Peyresq. We also thank Jose Luis Gruver for 
useful discussions, and Ana Maria Llois and Ruben Oscar Weht for reading the manuscript. 



[1] P. Ullersma, Physica 32, 27-55 (1966). 

[2] M.A. Castagnino, F.H. Gaioli, and E. Gunzig, Fund. Cosmic Phys. 16, 221-375 (1996) and 
references therein. 

[3] R. Davidson and J.J. Kozak, J. Math. Phys. 12, 903-917 (1971). 

[4] J.L. Gruver, J. Aliaga, H.A. Cerdeira, and A.N. Proto, Phys. Rev. E 51, 6263-6266 (1995). 

[5] K. Lindenberg and B.J. West, Phys. Rev. A 30, 568-582 (1984) and references therein. 

[6] K.O. Priedrichs, Commun. Pure Appl. Math. 1, 361-406 (1948). 

[7] F. Haake and R. Reibold, Phys. Rev. A 32, 2462-2475 (1985). 

[8] L. van Hove, Physica 21, 517-540 (1955). See also E.B. Davies, Commun. Math. Phys. 33, 
171-186 (1973). 

[9] E.C.G. Sudarshan, C.B. Chiu, and V. Gorini, Phys. Rev. D 18, 2914-2929 (1978). 

[10] M.A. Castagnino and R. Laura, Phys. Rev. A 56, 108-119 (1997). 

[11] P. Exner, Open Quantum Systems and Feynman Integrals (Reidel, Amsterdam, 1985). 

[12] L.A. Khalfin, Zh. Eksp. Teor. Fiz. 33, 1371-1382 (1957) [Sov. Phys. JETP 6, 1053-1063 (1958)]. 



40 



